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Abstract. Millisecond pulsars occur abundantly in globular clusters. They are expected to be 
responsible for several spectral components in the radio through 7 -ray waveband (e.g., involving 
synchrotron and inverse Compton emission), as have been seen by Radio Telescope Effelsberg, 
Chandra X-ray Observatory, Fermi Large Area Telescope, and the High Energy Stereoscopic 
System (H.E.S.S.) in the case of Terzan 5 (with fewer spectral components seen for other globular 
clusters). H.E.S.S. has recently performed a stacking analysis involving 15 non-detected globular 
clusters and obtained quite constraining average flux upper limits above 230 GeV. We present 
a model that assumes millisecond pulsars as sources of relativistic particles and predicts multi¬ 
wavelength emission from globular clusters. We apply this model to the population of clusters 
mentioned above to predict the average spectrum and compare this to the H.E.S.S. upper limits. 
Such comparison allows us to test whether the model is viable, leading to possible constraints 
on various average cluster parameters within this framework. 


1. Introduction 

There are nearly 160 Galactic globular clusters (GCs) known |14] . Each consists of hundreds of 
thousands of stars held together by their mutual gravity, and they orbit the Galactic Centre in a 
spherical distribution. GCs are typically about ten gigayears old and are therefore expected to 
harbour many evolved stellar objects, since the latter should have had ample time to complete 
their evolutionary processes. The high stellar densities in the cores of GCs also enhance stellar 
encounter rates, facilitating the formation of objects such as low-mass X-ray binaries (LMXRBs), 
cataclysmic variables, white dwarfs, and pulsars [T8j. LMXRBs are believed to be the progenitors 
of millisecond pulsars (MSPs; [5]), and since they occur abundantly in GCs, the same should 
hold true for MSPs. This is indeed found to be the case: 28 of the Galactic GCs contain more 
than 144 confirmed radio pulsar^, the bulk of these being MSPs. It was furthermore estimated 
that there should be 2 600 — 4 700 Galactic GC MSPs observable in 7 -rays [T]. 

GCs are multi-wavelength objects, being visible from radio to the highest energies. For 
example, several radio structures are visible in the direction of Terzan 5 at 11 cm and 21 cm [2j, 
while diffuse X-ray emission has also been detected from this GC lain]- The Fermi Large 
Area Telescope (LAT) plausibly detected about a dozen GCs in the GeV energy band [IT], and 
their spectral characteristics point to the cumulative emission from a population of GC MSPs. 


^ http://www.naic.edu/~pfreire/GCpsr.html 


In the TeV domain, H.E.S.S. has published upper limits for 47 Tucanae [3], and detected a 
very-high-energy (VHE) excess in the direction of Terzan 5 [2]. 

Several models have been proposed to explain the observed GC spectra. The total GeV 
contribution from GG MSPs was estimated by summing up individual predicted pulsed curvature 
radiation (CR) spectra from a population of MSPs [T3l[T9]. An alternative scenario [ 8 ] assumed 
that the GeV emission was due to inverse Compton (IC) radiation by leptons escaping from 
the MSP magnetospheres, upscattering Cosmic Microwave Background (CMB), stellar, as well 
as Galactic background (infrared and optical) photons. Another model [B] considered MSPs 
that accelerate leptons either at the shocks originating during collisions of neighbouring pulsar 
winds or inside the pulsar magnetospheres. These leptons escape from the magnetospheres and 
diffuse through the GC, encountering target photon fields such as optical and CMB emission. 
The latter are upscattered via the IC process, leading to GeV / TeV emission. This model was 
extended [ 201 123] using alternative particle injection spectra, and also calculating the expected 
synchrotron radiation (SR) from leptons moving in a homogeneous cluster magnetic held B. 
Recently, the model was signihcantly rehned |15] and now includes a line-of-sight calculation 
of the X-ray surface brightness (which is used to constrain the diffusion coefficient) as well as 
full particle transport, assuming spherical symmetry and a steady-state regime. Apart from 
leptonic models, there two other ideas. Non-accreting white dwarfs may make a considerable 
contribution to the cumulative 7 -ray hux seen from GGs |7], depending on their abundance. 
Lastly, a model invoking 7 -ray burst remnants as sources of energetic leptons and hadrons was 
also put forward [ 10 ] . A short burst may accelerate hadrons, which may in turn collide with 
ambient target nuclei, leading to particles that eventually decay into 7 -rays. In this model. 
X-rays may result from IG on optical stellar photons by primary electrons accelerated by the 
relativistic blast wave. 

We note that a search for diffuse X-ray emission from several Fermi-detected GGs [12] failed 
to detect any such emission above the Galactic background level, although a recent reanalysis 
of the archival Chandra data detected a new diffuse X-ray emission feature within the half-mass 
radius of 47 Tucanae [22]. In the context of the leptonic MSP models, these upper limits and 
detection place constraints on parameters such as the number of embedded MSPs Vmsp as well 
as the typical cluster held strength B. One way to help discriminate between the various models 
would be to follow a population approach. H.E.S.S. has recently obtained an upper limit to 
the average TeV hux of 15 non-detected GGs [3]. They noted that their upper limit was lower 
than the hux predicted by a simple leptonic scaling model by a factor of ~ 3 — 30, depending on 
model assumptions. This paper represents a hrst attempt to model the average TeV hux from 
the population of 15 GGs using our rehned leptonic GG model. Our motivation is to compare 
our results with those of the scaling model, and second, to assess the plausibility of the MSPs 
scenario by testing whether our average spectrum satishes the TeV upper limits. In Section |21 
we briehy discuss our model, after which follow a description of the assumed parameters and 
calculation method (Section |2D, as well as our conclusions (Section Hj). 

2. The basic GC model 

We solve the following transport equation numerically 

^ = v.(K.v„4-^(7~e)+e, (1) 

with Ug the electror@ density, which is a function of the radius vector 77 , Eg the electron energy, 
K, the diffusion tensor, Eg > 0 the particle energy losses, and Q the source term. We assume a 

^ We use the word ‘electron’ in a wider sense. Pair production may take place, additionally providing the 
possibility of positrons as efficient radiators. 



source term located at rg = that is of the form 


Q = Qo-Ee ^S{fs - fc), 


( 2 ) 


with the spectral index F. We assume stationarity and spherical symmetry, leading to 



( 3 ) 


with K the scalar diffusion coefficient. This may be rewritten as 


^ = —A (r -n^ + o] 

dEe Ee \ drs) ®(9£'e J 


( 4 ) 


and is solved numerically using a Crank-Nicolson algorithm. Our particle injection spectrum is 
normalised as follows: 

/* ^e.max 


EeQ dEg — .iW/[SP^(4^sd)) 


( 5 ) 


with rj the particle conversion efficiency, and {L^d) the average MSP spin-down luminosity. After 
solving for Ue, we calculate the SR and IC fluxes as detailed in m- We assume Bohm diffusion 


K{rs,Ee) 


cEe 

3eB{rs) ’ 


( 6 ) 


where c and e denote speed of light and electron charge. As an alternative, we also investigated 
a coefficient of the form ^ 

K{rs, Ee) = Ko{rs) , (7) 

with Egfi = 1 TeV and a = 0.6. For simplicity, we only considered spatially constant kq and B. 


3. Estimating the average GC spectrum - a first approach 

The model described in Section [2] has been successfully applied to the case of Terzan 5 [15] , 
where we performed a line-of-sight integration of the X-ray flux in order to constrain the lepton 
diffusion coefficient. We found that values of kq ~ 10“^ kpc^Myr”^ gave good hts. In this paper, 
we apply the same model to 15 non-detected GCs [3], using fixed parameters as noted in tabled] 
We have used the values of [T] for Nm^ where possible, and iV* values from m, and obtained 
distances d and structural parameterqf] Tc (core radius), rhm (half-light radius, used as a proxy 
for half-mass radius), and rt (tidal radius) from [T3]. Typical values of 77 ~ 0.01, Xmsp ~ 25, 
and (Lsd) ~ 2 x 10^^ ergs“^ led to values for the source strength Qo- We have used B = 5 fiG 
and T = 2.0 throughout, and used optical GC plus CMB photons as background fields for our 
IC calculation (assuming an average stellar temperature of T = 4 500 K). 

We hrst produced SR and IC spectra for each individual cluster. We next added these spectra 
and divided by the number of GCs to obtain the average predicted spectrum. This is indicated 
in hgure[T](differential flux) and figure [2](integral flux). None of the single-cluster spectra violate 
the TeV upper limits. The stacked upper limits are F{> 0.23 TeV) = 3.3 x 10“^^ cm“^s“^ for 
a point-like source analysis, and F{> 0.23 TeV) = 4.5 x 10“^^ cm“^s“^ for an extended source 
analysis |3|. Our average spectrum gives F{> 0.23 TeV) ~ 2 x 10“^^ cm“^s“^, which satishes 

® See also littp://gclusters.altervista.org/index.php 








Figure 1. Predicted differential spectra 
E‘^dN/dE for 15 single GCs (thin lines), 
as well as the average spectrum (thick red 
line). The two components represent the 
SR and IC spectra. 



Figure 2. Same as figurelH but for integral 
flux E{>E). Upper limits for the point-like 
source analysis are indicated by diamonds, 
and those for the extended source analysis 
by squares [3]. Associated spectra and 
upper limits have the same colour. The 
larger (lowest) symbols are for the average 
spectrum. 


these observational upper limits. However, we note that different choices of parameters will lead 
to average fluxes that could even exceed these limits in some cases. 

We next tested the effect of changing the diffusion coefficient. Figure [3] indicates the average 
spectra assuming Bohm diffusion (solid red line), a coefficient n = kq{E^/1 TeV)*^'® with 
Ko = 10““^ kpc^Myr“^ (dashed green line), and kq = 10“^ kpc^Myr“^ (dotted blue line). The 
change in spectral shape indicates the different energy dependencies {k oc E^'^ for Bohm diffusion 
versus k oc Fg ®). The corresponding integral fluxes are shown in figure 01 


Table 1. Assumed GC parameters. 


GG Name 

d 

(kpc) 

w 

(105) 

Amsp 

Qo 

(lO^^/erg/s) (') 

^hm 

(') 

n 

(') 

NGG 104 (47 Tuc) 

4.5 

4.57 

00+15 

9.55 

0.36 

3.17 

42.86 

NGG 6388 

9.9 

5.8 


52.1 

0.12 

0.52 

6.21 

NGG 7078 

10.4 

4.13 

25 (< 56) 

7.24 

0.14 

1.00 

21.5 

Terzan 6 

6.8 

0.29 

25 

7.24 

0.05 

0.44 

17.39 

Terzan 10 

5.8 

0.38 

25 

7.24 

0.9 

1.55 

5.06 

NGG 6715 

26.5 

4.79 

25 

7.24 

0.09 

0.82 

7.47 

NGG 362 

8.6 

1.58 

25 

7.24 

0.18 

0.82 

16.11 

Pal 6 

5.8 

0.31 

25 

7.24 

0.66 

1.2 

8.36 

NGG 6256 

10.3 

0.21 

25 

7.24 

0.02 

0.86 

7.59 

Djorg 2 

6.3 

1.0 

25 

7.24 

0.33 

1.0 

10.53 

NGG 6749 

7.9 

0.24 

25 

7.24 

0.62 

1.1 

5.21 

NGG 6144 

8.9 

0.48 

25 

7.24 

0.94 

1.63 

33.25 

NGG 288 

8.9 

0.32 

25 

7.24 

1.35 

2.23 

12.94 

HP 1 

8.2 

0.48 

25 

7.24 

0.03 

3.1 

8.22 

Terzan 9 

7.1 

0.02 

25 

7.24 

0.03 

0.78 

8.22 








Figure 3. Predicted average differential 
spectra. The line types indicate differ¬ 
ent diffusion coefficient assumptions as 
noted in the legend. 



Figure 4. Same as figure El but for 
the integrated flux. The red square 
and diamond indicate the extended and 
point-like H.E.S.S. upper limits. 


Our next task is a rigorous assessment of the error on the average integral spectrum, taking 
into account the range of values that each free parameter may assume, for each individual cluster. 
Preliminary investigations have shown that we can reduce the number of free parameters to five. 
Varying these parameters on a grid and calculating the average flux value and spread give 
quite large errors [ 2 T]. This process is very time-consuming, given the enormous number of 
combinations in which a cumulative (or average) spectrum may be obtained (15 GCs, with at 
least 5 free parameters each). 

4. Conclusions 

We have described the application of our leptonic GC model to a population of clusters that 
have been observed, but not detected, in VHE 7 -rays. The fact that upper limits are available 
presented us with the opportunity to test our model using a sizable sample of clusters. The 
upper limit for the stacked (average) flux is stringent, and a simple scaling model violated this 
limit. In contrast, none of the our single-cluster spectral predictions violated the respective 
upper limits. Our average model spectrum also satisfied the stacked upper limits for the given 
choice of parameters (tabled]), which we regarded as reasonable. There is, however, considerable 
uncertainty in the single-GG parameters, so that the intrinsic error on the predicted average 
spectrum will be quite large. As an example, we showed that different assumptions for the 
diffusion coefficient lead to significant changes in spectral and flux. One should therefore attempt 
to reduce the number of free model parameters and also constrain their values so as to decrease 
the predicted flux error. This will allow one to more robustly test the viability of the MSP 
scenario for explaining the (non-)detection of TeV flux from GGs. 
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